step0：Create 4 network

require(intergraph)

## Loading required package: intergraph

require(igraph)

## Loading required package: igraph

## 
## Attaching package: 'igraph'

## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum

## The following object is masked from 'package:base':
## 
##     union

require(RColorBrewer)

## Loading required package: RColorBrewer

require(gplots)

## Loading required package: gplots

## 
## Attaching package: 'gplots'

## The following object is masked from 'package:stats':
## 
##     lowess

require(ggsci)

## Loading required package: ggsci

library(RColorBrewer)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:igraph':
## 
##     as_data_frame, groups, union

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ape)

## 
## Attaching package: 'ape'

## The following objects are masked from 'package:igraph':
## 
##     degree, edges, mst, ring

library(sna)

## Loading required package: statnet.common

## 
## Attaching package: 'statnet.common'

## The following objects are masked from 'package:base':
## 
##     attr, order

## Loading required package: network

## 
## 'network' 1.18.0 (2022-10-05), part of the Statnet Project
## * 'news(package="network")' for changes since last version
## * 'citation("network")' for citation information
## * 'https://statnet.org' for help, support, and other information

## 
## Attaching package: 'network'

## The following objects are masked from 'package:igraph':
## 
##     %c%, %s%, add.edges, add.vertices, delete.edges, delete.vertices,
##     get.edge.attribute, get.edges, get.vertex.attribute, is.bipartite,
##     is.directed, list.edge.attributes, list.vertex.attributes,
##     set.edge.attribute, set.vertex.attribute

## sna: Tools for Social Network Analysis
## Version 2.7 created on 2022-05-09.
## copyright (c) 2005, Carter T. Butts, University of California-Irvine
##  For citation information, type citation("sna").
##  Type help(package="sna") to get started.

## 
## Attaching package: 'sna'

## The following objects are masked from 'package:ape':
## 
##     consensus, degree

## The following objects are masked from 'package:igraph':
## 
##     betweenness, bonpow, closeness, components, degree, dyad.census,
##     evcent, hierarchy, is.connected, neighborhood, triad.census

library(maps)

#.rs.restartR()

## Read advice.csv
country <- read.csv("country_metadata.csv", header = TRUE, as.is = TRUE, stringsAsFactors = FALSE)
goods <- read.csv("exp_goods_2009.csv", header = TRUE, as.is = TRUE, stringsAsFactors = FALSE)
services <- read.csv("exp_services_2009.csv", header = TRUE, as.is = TRUE, stringsAsFactors = FALSE)

# Creat the igraph object 
goods_net <- graph_from_edgelist(as.matrix(goods[,c(1,2)]), directed = TRUE)
services_net <- graph_from_edgelist(as.matrix(services[,c(1,2)]), directed = TRUE)

## Assign "USD" attribute to the edges
E(goods_net)$USD <- goods$USD
E(services_net)$USD <- services$USD

goods_top <- igraph::delete.edges(goods_net, which(E(goods_net)$USD < quantile(E(goods_net)$USD, prob = 0.75, na.rm = TRUE)))
services_top <- igraph::delete.edges(services_net, which(E(services_net)$USD < quantile(E(services_net)$USD, prob = 0.75, na.rm = TRUE)))


# add cols
attify <- function(n_net){
    igraph::V(n_net)$ISO.alpha3.code <- country$ISO.alpha3.code
    igraph::V(n_net)$Country <- country$Country
    igraph::V(n_net)$Continent <- country$Continent
    igraph::V(n_net)$GDP.per.capita <- country$GDP.per.capita
    igraph::V(n_net)$Region.1 <- country$Region.1
    igraph::V(n_net)$Region.2 <- country$Region.2
    igraph::V(n_net)$Population2010.OECD.estimate <- country$Population2010.OECD.estimate
    igraph::V(n_net)$Area.sqkm <- country$Area.sqkm
    igraph::V(n_net)$centroid.lon <- country$centroid.lon
    igraph::V(n_net)$centroid.lat <- country$centroid.lat
    igraph::V(n_net)$capital <- country$capital
    igraph::V(n_net)$capital.lon <- country$capital.lon
    igraph::V(n_net)$capital.lat <- country$capital.lat
  return(n_net)
}

goods_net2 <- attify(goods_net)
services_net2 <- attify(services_net)
goods_top2 <- attify(goods_top)
services_top2 <- attify(services_top)

# Create a color palette

col4 <- colorRampPalette(brewer.pal(8, 'Set3'))(length(unique(V(goods_net2)$Continent)))
continent_colors <- setNames(col4, unique(V(goods_net2)$Continent))

# Plotting the network
p_goods <- plot(goods_net2, 
                vertex.size = 0.002 * (strength(goods_net2, mode = "in", weights = log2(E(goods_net2)$USD) + 1)), 
                vertex.color = continent_colors[V(goods_net2)$Continent], 
                vertex.label = NA,
                edge.arrow.size = 0.1, 
                edge.curved = TRUE, 
                layout = as.matrix(cbind(country$centroid.lon, country$centroid.lat)),
                edge.width = 0.002 * log2(E(goods_net2)$USD),
                main = "Trade In Goods Between The Countries")

legend(-2.5, 1,
       legend = names(continent_colors),
       pch = 19,
       col = continent_colors,
       cex = 0.5)

# Polt for all service
col4 <- colorRampPalette(brewer.pal(8, 'Set3'))(length(unique(V(services_net2)$Continent)))
continent_colors <- setNames(col4, unique(V(services_net2)$Continent))

p_service <- plot(services_net2, 
                  vertex.size = 0.002*(strength(services_net2, mode = "in", weights = log2(E(services_net2)$USD))), 
                  vertex.color = continent_colors[V(services_net2)$Continent],
                  vertex.label = NA,
                  # vertex.label.cex = 0.5,
                  edge.arrow.size = 0.1, 
                  edge.curved = TRUE, 
                  layout = as.matrix(cbind(country$centroid.lon, country$centroid.lat)),
                  edge.width = 0.002*log2(E(services_net2)$USD),
                  main = "Trade In Services Between The Countries")

legend(-2.5, 1,
       legend = names(continent_colors),
       pch = 19,
       col = continent_colors,
       cex = 0.5)

# Polt for top goods
col4 <- colorRampPalette(brewer.pal(8, 'Set3'))(length(unique(V(goods_top2)$Continent)))
continent_colors <- setNames(col4, unique(V(goods_top2)$Continent))
p_goods_top <- plot(goods_top2, 
                    vertex.size = 0.004*(strength(goods_top2, mode = "in", weights = log2(E(goods_top2)$USD))),
                    vertex.color = continent_colors[V(goods_top2)$Continent],
                    vertex.label = NA,
                    edge.arrow.size = 0.1, 
                    edge.curved = TRUE, 
                    layout = as.matrix(cbind(country$centroid.lon, country$centroid.lat)),
                    edge.width = 0.008*log2(E(goods_top2)$USD),
                    main = "Goods Trading Relationships In Top 25%")

legend(-2.5, 1,
       legend = unique(V(goods_top2)$Continent),
       pch = 19,
       col = continent_colors,
       text.width = 1,
       cex = 0.5)

# Create a color palette
col4 <- colorRampPalette(brewer.pal(8, 'Set3'))(length(unique(V(goods_net2)$Continent)))
continent_colors <- setNames(col4, unique(V(goods_net2)$Continent))

# Plotting the network
plot(goods_net2, 
     vertex.size = 0.002 * (strength(goods_net2, mode = "in", weights = log2(E(goods_net2)$USD) + 1)), 
     vertex.color = continent_colors[V(goods_net2)$Continent],
     vertex.label = NA,
     edge.arrow.size = 0.1, 
     edge.curved = TRUE, 
     layout = as.matrix(cbind(country$centroid.lon, country$centroid.lat)),
     edge.width = 0.002 * log2(E(goods_net2)$USD),
     main = "Trade In Goods Between The Countries"
)

# Add a legend

legend(-2.5, 1,
       legend = unique(V(goods_net2)$Continent),
       pch = 19,
       col = continent_colors,
       text.width = 1,
       cex = 0.5)

# Polt for top service
col4 <- colorRampPalette(brewer.pal(8, 'Set3'))(length(unique(V(services_top2)$Continent)))
continent_colors <- setNames(col4, unique(V(services_top2)$Continent))


p_services_top <- plot(services_top2, 
                    vertex.size = 0.004*(strength(services_top2, mode = "in", weights = log2(E(services_top2)$USD))),
                    vertex.color = continent_colors[V(goods_net2)$Continent],
                    vertex.label = NA,
                    edge.arrow.size = 0.1, 
                    edge.curved = TRUE, 
                    layout = as.matrix(cbind(country$centroid.lon, country$centroid.lat)),
                    edge.width = 0.008*log2(E(services_top2)$USD),
                    main = "Sevice Trading Relationships In Top 25%")
legend(-2.5, 1,
       legend = unique(V(services_top2)$Continent),
       pch = 19,
       col = continent_colors,
       text.width = 1,
       cex = 0.5)

##1.Consider the overall metrics (average path length, transitivity, and reciprocity) for the top goods and top services networks. Compare these two networks to networks modelled on the originals, created with the Erdős–Rényi model and with the configuration model (i.e., you should create one E-R network and one configuration model network based on the top goods network, and one E-R network and one configuration model network based on the top services network). What do these comparisons tell you about the nature of trade flows between countries? How do the two empirical networks differ? In your answer, make sure to define each of the metrics and give an intuitive interpretation for them.

# Generate an 100 times random sample for goods comparison

goods_ran <- lapply(rep(1, 100), function(x)
  sample_gnp(n=vcount(goods_top2), p=graph.density(goods_top2), directed = TRUE))

goods_ran_apl <- sapply(goods_ran, average.path.length)
goods_ran_tra <- sapply(goods_ran, transitivity)
goods_ran_rec <- sapply(goods_ran, reciprocity)


# Create 100 networks using Configuration model 

ind <- igraph::degree(goods_top2, mode='in')
outd <- igraph::degree(goods_top2, mode='out')

goods_con <- lapply(rep(1, 100), function(x)
  sample_degseq(out.deg=outd, in.deg=ind, method="simple"))

goods_con_apl <- sapply(goods_con, average.path.length)
goods_con_tra <- sapply(goods_con, transitivity)
goods_con_rec <- sapply(goods_con, reciprocity)

goods_table <- data.frame(c('average.path.length', 'transitivity','reciprocity'),
                       c(average.path.length(goods_top2), transitivity(goods_top2), reciprocity(goods_top2)),
                       c(mean(goods_ran_apl), mean(goods_ran_tra), mean(goods_ran_rec)),
                       c(mean(goods_con_apl), mean(goods_con_tra), mean(goods_con_rec)))
colnames(goods_table) <- c('Name', 'Top Goods','Top Goods Random','Top Goods Constructed')
goods_table

##                  Name Top Goods Top Goods Random Top Goods Constructed
## 1 average.path.length 1.8456002        1.8026006             2.0129183
## 2        transitivity 0.5871980        0.3579854             0.5081025
## 3         reciprocity 0.7564153        0.2002331             0.3318143

# Generate an 100 times random sample  for service comparison
services_ran <- lapply(rep(1, 100), function(x)
  sample_gnp(n=vcount(services_top2), p=graph.density(services_top2), directed = TRUE))

services_ran_apl <- sapply(services_ran, average.path.length)
services_ran_tra <- sapply(services_ran, transitivity)
services_ran_rec <- sapply(services_ran, reciprocity)

# Create 100 networks using Configuration model 

ind <- igraph::degree(services_top2, mode='in')
outd <- igraph::degree(services_top2, mode='out')

services_con <- lapply(rep(1, 100), function(x)
  sample_degseq(out.deg=outd, in.deg=ind, method="simple"))

services_con_apl <- sapply(services_con, average.path.length)
services_con_tra <- sapply(services_con, transitivity)
services_con_rec <- sapply(services_con, reciprocity)


services_table <- data.frame(c('average.path.length', 'transitivity','reciprocity'),
                             c(average.path.length(services_top2), transitivity(services_top2), reciprocity(services_top2)),
                             c(mean(services_ran_apl), mean(services_ran_tra), mean(services_ran_rec)),
                             c(mean(services_con_apl), mean(services_con_tra), mean(services_con_rec)))

colnames(services_table) <- c('Name', 'Top Services','Top Services Random','Top Services Configuration')
services_table

##                  Name Top Services Top Services Random
## 1 average.path.length    1.7485525           1.7500857
## 2        transitivity    0.6051970           0.4374527
## 3         reciprocity    0.8379994           0.2496224
##   Top Services Configuration
## 1                  1.9505484
## 2                  0.5408152
## 3                  0.3609007

# Answer：
# What do these comparisons tell you about the nature of trade flows between countries? 

# How do the two empirical networks  differ? 

# In  your  answer,  make  sure  to  define  each  of  the  metrics  and  give  an  intuitive interpretation for them.

##2.Which do you see as the most influential countries in the trade networks? Identify two potential meanings of “influence,” as proxied by different centrality measures. Justify your choice of each centrality measure. In that justification, present clear interpretations of what each centrality measure is capturing about the position of the countries. Calculate both of your chosen centrality measures on the top goods and top services networks and identify the top ten countries that has the highest values for each. There are likely to be many countries repeatedly identified across both, so it may be informative to consider the countries that appear uniquely as having high centrality in one network but not the other. Discuss the patterns you see and what it implies about the top countries and about trade flows, generally. Make explicit reference in your response to the concepts of social capital and brokerage/structural holes covered in the course material. [Note that if you use edge weight for a centrality measure, you need to identify how the calculation interprets those weights (i.e., do higher values mean greater closeness or greater distance?); in those cases where the measure assumes that higher values means greater distance, you should use 1/weight in the calculation.]

cents <- function(net){
  
   # Store current scipen value
  old_scipen <- options("scipen")$scipen
  # Set scipen to a high value to prevent scientific notation
  options(scipen = 999)
  
  bet <- igraph::betweenness(net, directed = FALSE, weight = E(net)$USD) 
  clo <- igraph::closeness(net, mode = "all", weights = 1/E(net)$weight)
  all <- data.frame(bet,pr)
  return(all)
}


# top 10 betweenness and closeness for goods

clo <- igraph::closeness(goods_top2, mode = "all", weights = 1/E(goods_top2)$USD) %>% 
  as.data.frame() 
colnames(clo) <-"goods_closeness"
arrange(clo, desc(goods_closeness)) %>% head(10)

##     goods_closeness
## USA         6929349
## CHN         6921370
## DEU         6899532
## CAN         6897758
## MEX         6888568
## HKG         6878391
## JPN         6877432
## GBR         6849432
## KOR         6842573
## FRA         6839074

bet <- igraph::betweenness(goods_top2, directed = TRUE, weight = E(goods_top2)$USD) %>% 
  as.data.frame() 
colnames(bet) <-"goods_betweenness"
arrange(bet, desc(goods_betweenness)) %>% head(10)

##     goods_betweenness
## CHN              1298
## USA              1257
## DEU              1250
## IND               955
## ZAF               857
## BEL               822
## CAN               803
## NLD               766
## CHE               734
## GBR               694

# top 10 betweenness and closeness for services
clo <- igraph::closeness(services_top2, mode = "all", weights = 1/E(services_top2)$USD) %>% 
  as.data.frame() 
colnames(clo) <-"services_closeness"
arrange(clo, desc(services_closeness)) %>% head(10)

##     services_closeness
## USA            2121266
## GBR            2111208
## DEU            2108010
## CAN            2106747
## JPN            2104786
## BMU            2100090
## FRA            2099151
## IRL            2097767
## NLD            2088290
## CHN            2088001

bet <- igraph::betweenness(services_top2, directed = TRUE, weight = E(services_top2)$USD) %>% 
  as.data.frame() 
colnames(bet) <-"services_betweenness"
arrange(bet, desc(services_betweenness)) %>% head(10)

##     services_betweenness
## USA                 1295
## IND                 1035
## CHN                 1011
## BEL                  831
## KOR                  827
## SGP                  827
## FRA                  793
## JPN                  780
## ESP                  774
## DEU                  751

#3.How does geography influence trade flows? Calculate assortativity by continent and by GDP on the top goods and top services networks. Use the spinglass community detection algorithm to divide both the goods and services networks into communities (make sure to consider edge weight and use set.seed(123456) to ensure that we get equivalent results). Plot the goods and services networks with nodes positioned by their latitude/longitude (so that your network should look roughly like a map of the world), and nodes coloured by their community membership. Make sure that your plots are legible and informative, with a legend. Discuss what each of these results implies about trade flows.

# top goods betweeness

as_gC <- assortativity.nominal(goods_top2, as.factor(V(goods_top2)$Continent), directed = T)
as_gG <- assortativity(goods_top2, as.factor(V(goods_top2)$GDP.per.capita), directed = T)
as_sC <- assortativity.nominal(services_top2, as.factor(V(services_top2)$Continent),directed = T)
as_sG <- assortativity(services_top2, as.factor(V(services_top2)$GDP.per.capita), directed = T)

assortativity_table <- data.frame(c('By Continent','By GDP'),
                                  c(as_gC,as_gG),
                                  c(as_sC,as_sG))
colnames(assortativity_table) <- c('assortativity','Goods','Service')
assortativity_table

##   assortativity        Goods    Service
## 1  By Continent -0.001593921  0.0590188
## 2        By GDP  0.002879183 -0.1775439

set.seed(123456)
goods_component <- decompose.graph(goods_top2)[[1]]
sp_g <- cluster_spinglass(goods_component, weights = E(goods_component)$USD)
sp_g$csize

## [1] 55 12 46 43

split(sp_g$names, sp_g$membership)

## $`1`
##  [1] "ARM" "AUT" "AZE" "BEL" "BFA" "BGR" "BIH" "BLR" "CHE" "CIV" "CMR" "CZE"
## [13] "DEU" "DNK" "DZA" "EGY" "ESP" "EST" "FIN" "FRA" "GBR" "GEO" "GRC" "HRV"
## [25] "HUN" "IRL" "ITA" "JOR" "KAZ" "LBN" "LTU" "LUX" "LVA" "MKD" "MLT" "MNE"
## [37] "NLD" "NOR" "POL" "ROU" "RUS" "SVK" "SVN" "SWE" "TUN" "TUR" "UKR" "CPV"
## [49] "MRT" "PRT" "CYP" "ISL" "MAR" "MDA" "SEN"
## 
## $`2`
##  [1] "BEN" "BWA" "MOZ" "TGO" "ZAF" "ZWE" "GHA" "NAM" "ZMB" "MLI" "MWI" "LSO"
## 
## $`3`
##  [1] "AGO" "ALB" "ATG" "BOL" "BRA" "BRB" "CAN" "CHN" "COL" "CRI" "DOM" "GTM"
## [13] "HKG" "HND" "ISR" "KGZ" "KHM" "KNA" "MEX" "NER" "NGA" "NIC" "PER" "SLV"
## [25] "SWZ" "URY" "USA" "CHL" "COG" "ECU" "GMB" "MAC" "TTO" "VEN" "ARG" "BMU"
## [37] "PRY" "VCT" "BLZ" "GUY" "JAM" "LCA" "MNG" "PSE" "SLB" "SUR"
## 
## $`4`
##  [1] "AFG" "ARE" "AUS" "BHR" "BRN" "ETH" "IDN" "IND" "IRN" "IRQ" "JPN" "KEN"
## [13] "KOR" "LKA" "MDG" "MUS" "MYS" "NPL" "NZL" "PAK" "PHL" "SAU" "SGP" "THA"
## [25] "TWN" "TZA" "UGA" "YEM" "BGD" "BHS" "RWA" "BDI" "KWT" "OMN" "PAN" "QAT"
## [37] "FJI" "LAO" "MDV" "MMR" "NCL" "VNM" "WSM"

plot(goods_component, 
     vertex.size = 0.004*(strength(goods_component, mode = "in", weights = log2(E(goods_component)$USD))),
     vertex.color = membership(sp_g),
     edge.width = 0.005*log(E(goods_component)$USD),
     edge.curved = TRUE, 
     vertex.label = NA,
     edge.arrow.size = 0.1,
     layout = as.matrix(cbind(country$centroid.lon, country$centroid.lat)),
     main = "Goods Spinglass Communities")

## Warning in layout[, 1] + label.dist * cos(-label.degree) * (vertex.size + :
## longer object length is not a multiple of shorter object length

## Warning in layout[, 2] + label.dist * sin(-label.degree) * (vertex.size + :
## longer object length is not a multiple of shorter object length

legend(-2.5,
       1,
       legend = c(1:4),
       pch = 19,
       cex = 0.5,
       col = categorical_pal(8)[c(1:4)])

# For services network
set.seed(123456)
services_component <- igraph::decompose.graph(services_top2)[[1]]
sp_s <- cluster_spinglass(services_component, weights = E(services_component)$USD)
sp_s$csize

## [1] 47 51 59

split(sp_s$names, sp_s$membership)

## $`1`
##  [1] "AGO" "ARG" "ARM" "ATG" "BDI" "BEN" "BHS" "BLZ" "BMU" "BOL" "BRA" "BRB"
## [13] "BWA" "CAF" "CAN" "CHL" "COG" "COL" "CRI" "DOM" "ECU" "FJI" "GEO" "GTM"
## [25] "GUY" "HND" "ISR" "JAM" "KNA" "LCA" "LSO" "MEX" "MWI" "NIC" "PAN" "PER"
## [37] "PRY" "PSE" "RWA" "SLV" "TTO" "URY" "USA" "VCT" "VEN" "ZMB" "ZWE"
## 
## $`2`
##  [1] "AFG" "ARE" "AUS" "AZE" "BGD" "BHR" "BRN" "CHN" "CMR" "EGY" "ETH" "HKG"
## [13] "IDN" "IND" "IRN" "IRQ" "JOR" "JPN" "KEN" "KGZ" "KHM" "KOR" "KWT" "LAO"
## [25] "LBN" "LKA" "MAC" "MDG" "MDV" "MMR" "MNG" "MOZ" "MUS" "MYS" "NCL" "NPL"
## [37] "NZL" "OMN" "PAK" "PHL" "QAT" "SAU" "SGP" "THA" "TON" "TWN" "TZA" "UGA"
## [49] "VNM" "WSM" "YEM"
## 
## $`3`
##  [1] "ALB" "AUT" "BEL" "BFA" "BGR" "BIH" "BLR" "CHE" "CIV" "CPV" "CYP" "CZE"
## [13] "DEU" "DNK" "DZA" "ESP" "EST" "FIN" "FRA" "GBR" "GHA" "GMB" "GRC" "HRV"
## [25] "HUN" "IRL" "ISL" "ITA" "KAZ" "LTU" "LUX" "LVA" "MAR" "MDA" "MKD" "MLI"
## [37] "MLT" "MNE" "MRT" "NAM" "NER" "NGA" "NLD" "NOR" "POL" "PRT" "ROU" "RUS"
## [49] "SEN" "SUR" "SVK" "SVN" "SWE" "SWZ" "TGO" "TUN" "TUR" "UKR" "ZAF"

plot(services_component, 
     vertex.size = 0.004*(strength(services_component, mode = "in", weights = log2(E(services_component)$USD))),
     vertex.color = membership(sp_s),
     edge.width = 0.005*log(E(services_component)$USD),
     edge.curved = TRUE, 
     vertex.label = NA,
     edge.arrow.size = 0.1,
     layout = as.matrix(cbind(country$centroid.lon, country$centroid.lat)),
     main = "Services Spinglass Communities")

## Warning in layout[, 1] + label.dist * cos(-label.degree) * (vertex.size + :
## longer object length is not a multiple of shorter object length

## Warning in layout[, 1] + label.dist * cos(-label.degree) * (vertex.size + :
## longer object length is not a multiple of shorter object length

legend(-2.5,
       1,
       legend = c(1:3),
       pch = 19,
       cex = 0.5,
       col = categorical_pal(8)[c(1:3)])

4.Which countries seem to have similar approaches to trading in goods and services? Calculate the structural equivalence (separately) on the top goods and top services networks (using the equiv.clust() function). Plot the dendrograms that result from each. Plot the top goods and top services networks again with nodes positioned as above by their latitude/longitude, and now with nodes coloured by the resulting structural equivalence classes (where you should divide each into four groups). Make sure that your plots are legible and informative, with a legend. Describe and compare the resulting blocks (for this, you will want to look at the dendrogram, at the memberships, and potentially at block models or mixing matrices based on the equivalence classes) and discuss what each of these results imply about global trade relations.

# Calculate and plot the structural equivalence on top goods

sna_goods_top <-intergraph::asNetwork(goods_top2)
st_eq_goods <- equiv.clust(sna_goods_top,mode = "digraph")
plot(st_eq_goods, cex = 0.24, main = "Top Goods Cluster")
rect.hclust(st_eq_goods$cluster, k = 4)

bm_goods <- blockmodel(sna_goods_top, ec = st_eq_goods, k = 4)
bm_goods

## 
## Network Blockmodel:
## 
## Block membership:
## 
## AFG AGO ALB ARE ARM ATG AUS AUT AZE BEL BEN BFA BGR BHR BIH BLR BOL BRA BRB BRN 
##   1   1   1   2   1   1   2   2   1   3   1   1   1   4   1   1   1   3   1   1 
## BWA CAN CHE CHN CIV CMR COL CRI CZE DEU DNK DOM DZA EGY ESP EST ETH FIN FRA GBR 
##   1   3   3   3   1   1   4   1   2   3   2   1   4   4   3   1   1   2   3   3 
## GEO GRC GTM HKG HND HRV HUN IDN IND IRL IRN IRQ ISR ITA JOR JPN KAZ KEN KGZ KHM 
##   1   2   1   2   1   1   2   2   3   2   4   4   2   3   4   3   1   4   1   1 
## KNA KOR LBN LKA LTU LUX LVA MDG MEX MKD MLT MNE MOZ MUS MYS NER NGA NIC NLD NOR 
##   1   3   4   1   1   1   1   1   4   1   1   1   1   1   2   1   4   1   3   2 
## NPL NZL PAK PER PHL POL ROU RUS SAU SGP SLV SVK SVN SWE SWZ TGO THA TUN TUR TWN 
##   1   4   4   4   4   2   2   3   2   2   1   2   1   2   1   1   2   4   3   2 
## TZA UGA UKR URY USA YEM ZAF ZWE BGD BHS CHL COG CPV ECU GHA GMB MAC MRT NAM PRT 
##   4   1   2   1   3   4   2   1   4   1   4   1   1   4   1   1   1   1   1   2 
## RWA STP TTO VEN ZMB ARG BDI BMU CYP ISL KWT MAR MDA MLI MWI OMN PAN PRY QAT VCT 
##   1   1   1   4   1   4   1   1   1   1   4   4   1   1   1   4   4   1   4   1 
## BLZ CAF FJI GUY JAM KIR LAO LCA LSO MDV MMR MNG NCL PSE SEN SLB SUR TON VNM WSM 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   2   1 
## 
## Reduced form blockmodel:
## 
##   AFG AGO ALB ARE ARM ATG AUS AUT AZE BEL BEN BFA BGR BHR BIH BLR BOL BRA BRB BRN BWA CAN CHE CHN CIV CMR COL CRI CZE DEU DNK DOM DZA EGY ESP EST ETH FIN FRA GBR GEO GRC GTM HKG HND HRV HUN IDN IND IRL IRN IRQ ISR ITA JOR JPN KAZ KEN KGZ KHM KNA KOR LBN LKA LTU LUX LVA MDG MEX MKD MLT MNE MOZ MUS MYS NER NGA NIC NLD NOR NPL NZL PAK PER PHL POL ROU RUS SAU SGP SLV SVK SVN SWE SWZ TGO THA TUN TUR TWN TZA UGA UKR URY USA YEM ZAF ZWE BGD BHS CHL COG CPV ECU GHA GMB MAC MRT NAM PRT RWA STP TTO VEN ZMB ARG BDI BMU CYP ISL KWT MAR MDA MLI MWI OMN PAN PRY QAT VCT BLZ CAF FJI GUY JAM KIR LAO LCA LSO MDV MMR MNG NCL PSE SEN SLB SUR TON VNM WSM 
##            Block 1    Block 2   Block 3   Block 4
## Block 1 0.01672625 0.07605877 0.2894911 0.0276886
## Block 2 0.11884183 0.85384615 0.9796380 0.5123626
## Block 3 0.38797092 0.99547511 1.0000000 0.9432773
## Block 4 0.03892456 0.31730769 0.7521008 0.2407407

# plot network
plot(goods_top2, 
     vertex.size = 0.004*(strength(goods_top2, mode = "in", weights = log2(E(goods_top2)$USD))),
     vertex.color = bm_goods$block.membership[order(bm_goods$order.vector)],
     edge.width = 0.005*log(E(goods_top2)$USD),
     edge.curved = TRUE, 
     vertex.label = NA,
     edge.arrow.size = 0.1,
     layout = as.matrix(cbind(country$centroid.lon, country$centroid.lat)),
     main = "Top Goods Structural Equivalence")
legend(-2.5,
       1,
       legend = c(1:4),
       pch = 19,
       cex = 0.5,
       col = categorical_pal(8)[c(1:4)])

# Calculate and plot the structural equivalence on top services
sna_services_top <- intergraph::asNetwork(services_top2)
st_eq_services <- equiv.clust(sna_services_top,mode = "digraph")
plot(st_eq_services, cex = 0.24, main = "Top Services Cluster")
rect.hclust(st_eq_services$cluster, k = 4)

bm_services <- blockmodel(sna_services_top, ec = st_eq_services, k = 4)

# for clear membership
tp_1 <- as.data.frame(bm_services$block.membership[order(bm_services$order.vector)])
service_memb_list <- cbind(country$Country, tp_1)
colnames(service_memb_list) <- c("country","cluster")
arrange(service_memb_list,cluster)

##                              country cluster
## 1                        Afghanistan       1
## 2                            Albania       1
## 3                            Armenia       1
## 4                Antigua and Barbuda       1
## 5                            Burundi       1
## 6                              Benin       1
## 7                       Burkina Faso       1
## 8                         Bangladesh       1
## 9                            Bahamas       1
## 10            Bosnia and Herzegovina       1
## 11                            Belize       1
## 12                           Bermuda       1
## 13                           Bolivia       1
## 14                          Barbados       1
## 15                 Brunei Darussalam       1
## 16                          Botswana       1
## 17          Central African Republic       1
## 18                     Cote d'Ivoire       1
## 19                          Cameroon       1
## 20                             Congo       1
## 21                        Cabo Verde       1
## 22                        Costa Rica       1
## 23                    Dominican Rep.       1
## 24                           Ecuador       1
## 25                          Ethiopia       1
## 26                              Fiji       1
## 27                           Georgia       1
## 28                             Ghana       1
## 29                            Gambia       1
## 30                         Guatemala       1
## 31                            Guyana       1
## 32                          Honduras       1
## 33                           Iceland       1
## 34                           Jamaica       1
## 35                             Kenya       1
## 36                        Kyrgyzstan       1
## 37                          Cambodia       1
## 38                          Kiribati       1
## 39             Saint Kitts and Nevis       1
## 40            Lao People's Dem. Rep.       1
## 41                       Saint Lucia       1
## 42                         Sri Lanka       1
## 43                           Lesotho       1
## 44               Republic of Moldova       1
## 45                        Madagascar       1
## 46                          Maldives       1
## 47                   North Macedonia       1
## 48                              Mali       1
## 49                           Myanmar       1
## 50                        Montenegro       1
## 51                          Mongolia       1
## 52                        Mozambique       1
## 53                        Mauritania       1
## 54                         Mauritius       1
## 55                            Malawi       1
## 56                           Namibia       1
## 57                     New Caledonia       1
## 58                             Niger       1
## 59                         Nicaragua       1
## 60                             Nepal       1
## 61                          Paraguay       1
## 62                State of Palestine       1
## 63                            Rwanda       1
## 64                           Senegal       1
## 65                   Solomon Islands       1
## 66                       El Salvador       1
## 67             Sao Tome and Principe       1
## 68                          Suriname       1
## 69                          Eswatini       1
## 70                              Togo       1
## 71                             Tonga       1
## 72               Trinidad and Tobago       1
## 73                           Tunisia       1
## 74           United Rep. of Tanzania       1
## 75                            Uganda       1
## 76                           Uruguay       1
## 77  Saint Vincent and the Grenadines       1
## 78                             Samoa       1
## 79                             Yemen       1
## 80                            Zambia       1
## 81                          Zimbabwe       1
## 82                            Angola       2
## 83                         Argentina       2
## 84                           Austria       2
## 85                        Azerbaijan       2
## 86                          Bulgaria       2
## 87                           Bahrain       2
## 88                           Belarus       2
## 89                            Brazil       2
## 90                             Chile       2
## 91                          Colombia       2
## 92                            Cyprus       2
## 93                           Czechia       2
## 94                           Algeria       2
## 95                             Egypt       2
## 96                           Estonia       2
## 97                           Finland       2
## 98                            Greece       2
## 99                           Croatia       2
## 100                          Hungary       2
## 101                        Indonesia       2
## 102                          Ireland       2
## 103              Iran (Islamic Rep.)       2
## 104                             Iraq       2
## 105                           Israel       2
## 106                           Jordan       2
## 107                       Kazakhstan       2
## 108                           Kuwait       2
## 109                          Lebanon       2
## 110                        Lithuania       2
## 111                       Luxembourg       2
## 112                           Latvia       2
## 113                 China, Macao SAR       2
## 114                          Morocco       2
## 115                           Mexico       2
## 116                            Malta       2
## 117                         Malaysia       2
## 118                          Nigeria       2
## 119                           Norway       2
## 120                      New Zealand       2
## 121                             Oman       2
## 122                         Pakistan       2
## 123                           Panama       2
## 124                             Peru       2
## 125                      Philippines       2
## 126                           Poland       2
## 127                         Portugal       2
## 128                            Qatar       2
## 129                          Romania       2
## 130               Russian Federation       2
## 131                         Slovakia       2
## 132                         Slovenia       2
## 133                           Sweden       2
## 134                           Turkey       2
## 135                           Taiwan       2
## 136                          Ukraine       2
## 137       Venezuela (Boliv. Rep. of)       2
## 138                         Viet Nam       2
## 139                     South Africa       2
## 140             United Arab Emirates       3
## 141                        Australia       3
## 142                          Belgium       3
## 143                           Canada       3
## 144                      Switzerland       3
## 145                          Denmark       3
## 146                            Spain       3
## 147             China, Hong Kong SAR       3
## 148                            India       3
## 149                            Italy       3
## 150                            Japan       3
## 151                    Rep. of Korea       3
## 152                     Saudi Arabia       3
## 153                        Singapore       3
## 154                         Thailand       3
## 155                            China       4
## 156                          Germany       4
## 157                           France       4
## 158                   United Kingdom       4
## 159                      Netherlands       4
## 160                    United States       4

plot(services_top2, 
     vertex.size = 0.004*(strength(services_top2, mode = "in", weights = log2(E(services_top2)$USD))),
     vertex.color = bm_services$block.membership[order(bm_services$order.vector)],
     edge.width = 0.005*log(E(services_top2)$USD),
     edge.curved = TRUE, 
     vertex.label = NA,
     edge.arrow.size = 0.1,
     layout = as.matrix(cbind(country$centroid.lon, country$centroid.lat)),
     main = "Top Services Structural Equivalence")
legend(-2.5,
       1,
       legend = c(1:4),
       pch = 19,
       cex = 0.5,
       col = categorical_pal(8)[c(1:4)])

#5

last_one<- data.frame(c('log_dist','log_in_gdp','log_out_gdp','log_in_pop','log_out_pop','mutual'),
           exp(c(-1.2,0.83,1.2,0.82,1.08,1.49)),
           exp(c(-0.91,1.10,1.34,0.88,0.79,2.12)))
colnames(last_one) <- c("odds  ratios","goods","service")

last_one

##   odds  ratios     goods   service
## 1     log_dist 0.3011942 0.4025242
## 2   log_in_gdp 2.2933187 3.0041660
## 3  log_out_gdp 3.3201169 3.8190435
## 4   log_in_pop 2.2704998 2.4108997
## 5  log_out_pop 2.9446796 2.2033964
## 6       mutual 4.4370955 8.3311375